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We develop an analytic approach to calculation of the temperature and polarisa- 
■ tion power spectra of the cosmic microwave background due to inflationary gravi- 

tational waves. This approach complements the more precise numerical results by 



providing insight into the physical origins of the features in the power spectra. 
We explore the use of analytic approximations for the gravitational-wave evolution, 
making use of the WKB approach to handle the radiation-matter transition. In the 
0^ \ process, we describe scaling relations for the temperature and polarisation power 

\ spectra. We illustrate the dependence of the amplitude, shape, and peak locations 

on the details of recombination, the gravitational-wave power spectrum, and the 
^ , cosmological parameters, and explain the origin of the peak locations in the tem- 

\ perature and polarisation power spectra. The decline in power on small scales in the 

polarisation power spectra is discussed in terms of phase-damping. In an appendix 
we detail numerical techniques for integrating the gravitational-wave evolution in 
the presence of anisotropic stress from free-streaming neutrinos. 
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1 Introduction 



The evidence from the cosmic microwave background (CMB) for a flat Uni- 
verse [1,2,3,4,5,6,7,8,9,10,11,12] and a nearly scale-invariant spectrum of pri- 
mordial adiabatic perturbations [2,3,4,5,6,7,8,9,10,11,12], in good agreement 
with the predictions of inflation [13,14,15,16,17,18], motivates additional ob- 
servational probes of inflation. One such probe is the polarisation signature 
[19,20,21,22] of the stochastic gravitational- wave background [23,24,25,26] pro- 
duced during inflation, which has now become the target of several ground- 
based experiments [27,28,29,30], as well as an Einstein vision experiment in 
NASA's science roadmap [31]. 

Large-angle CMB temperature fluctuations from these gravitational waves 
(tensor metric perturbations) were flrst considered in Refs. [23,24,25,26] and 
the polarisation was flrst considered in Ref. [32]. Now, the most precise predic- 
tions for these power spectra come from numerical calculations [33]. Like the 
power spectra for density perturbations (scalar metric perturbations), which 
exhibit wiggles due to acoustic waves in the primordial baryon-photon fluid, 
the temperature and polarisation power spectra from gravitational waves ex- 
hibit wiggles due to oscillations of tensor modes as they enter the horizon. 
The wiggles in the density-perturbation power spectra were predicted origi- 
nally by Sunyaev and Zeldovich [34] and Peebles and Yu [35], and explained 
later elegantly with a semi-analytic approach in a paper by Hu and Sugiyama 
[36]. 

The goal of this paper is to present an analytic account of the features in 
the tensor power spectra. Such an approach explains the origin of the fea- 
tures in the temperature and polarisation power spectra and illustrates the 
dependence of these features on the tensor power spectrum, cosmological pa- 
rameters, and details of the recombination history. The intuition provided by 
such an approach complements the more precise results of numerical calcula- 
tions. In particular, we explain here the location of the wiggles in the tensor 
temperature and polarisation power spectra, and why the bumps in the curl 
component of the polarisation are smoother than those in the curl-free com- 
ponent. We also show how the amplitude of the polarisation depends on the 
details of the recombination history. Our approach is analogous to that for 
scalar modes given in Ref. [36]. We discuss how measurement of the locations 
of these peaks can provide an independent probe of cosmological parameters. 

The organisation of the paper is as follows. In Section 2, we write the exact 
equations for CMB fluctuations from tensor perturbations. The exact equa- 
tions consist of the evolution of the gravitational waves, the visibility function, 
the source function, and projection factors. Next we develop a qualitative 
understanding of the physics contained in these relations in Section 3. The 
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remainder of the paper then investigates individually each ingredient in the 
exact calculation. Section 4 discusses the evolution of the gravitational-wave 
perturbation. Section 5 discusses the effect of the recombination history on 
the power spectrum. Section 6 discusses the projection factors, and Section 7 
the source function. Finally we comment on the dependence on cosmological 
parameters and detectability in Section 8. We include two Appendices that 
discuss the numerical techniques required to evolve the gravitational-wave am- 
plitude in the presence of neutrino anisotropic stress (Appendix A), and the 
application of the WKB approach to gravitational waves evolving through the 
matter-radiation transition (Appendix B). 



2 Exact Equations 

Here we present the exact equations required to evaluate the CMB power 
spectra from gravitational waves [19,21,37]. For simplicity, we will restrict 
ourselves to the case of a fiat FRW universe. Our emphasis will be on small 
scale structure and so reionisation and its effects on large scales will not be 
discussed; see Refs. [38,39] for more details on this topic. 

To provide the framework for temperature and polarisation anisotropics, we 
follow the formalism of Ref. [21]. For two other useful introductions into the 
subject see Refs. [40,41]. The CMB radiation field is characterised by the 
Stokes parameters I,Q, and U. The fourth Stokes parameter, V, is not gener- 
ated by Thomson scattering, and while it can be generated after last scattering, 
the expected amplitudes are small [42] and so can be neglected. 

While convenient, the Stokes parameters Q and U describing polarisation 
suffer from being co-ordinate dependent. Under a right-handed rotation by an 
angle ip in the plane perpendicular to the direction h of propagation, Q and 
U transform according to 

g' = (5cos2V' + C/sin2V', (1) 

[/' = -gsin2V' + C/cos2z^, (2) 
where e'j^ = cos'^ei -|- sin ■062 and e'j = — sin^/'ei + cos ^§2. These transfor- 
mation laws motivate the combinations Q ± ill which have definite spin-2. 
We can then expand these quantities in terms of the spin-weighted spherical 
harmonics which are appropriate to describe the projection of these quantities 
onto the unit sphere, 

T{n)^J2<'T,imYim{n), (3) 

Im 

{Q ± iU){h) = a±2M ±2Yim{fi), (4) 

Im 
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where Yijn{n) and ±2Yim{'n) are the spin-0 and spin-2 spin-weighted spher- 
ical harmonics. These expressions may be inverted to obtain the spherical- 
harmonic expansion coefficients, 

aTM- I dnY;jn)T{h), (5) 

a±2M = I dn±2Y;^ih){Q±iU){n). (6) 

Rather than work in terms of a±2,im it is advantageous to define two rotation- 
ally invariant quantities E and B by the relations 

fl£,im — —{0'2,lm + fl-2,/m)/2, (7) 



aB,lm — i{0'2,lm — a-2,im)/2. (8) 

These two quantities are equivalent to the curl and grad modes defined in Ref. 
[19]. The E mode is invariant under the parity transformation, while the B 
mode transforms with odd parity. 

Prom the above ax,im, we can form a series of correlation functions that char- 
acterise the statistics of the CMB perturbations. Of the six possible combi- 
nations, the TB and EE cross-correlations will vanish unless parity is some- 
how violated in the early Universe. The power spectra are defined as the 
rotationally-invariant quantities, 

Cxx'i = 21 + 1 '^^^xm'^^'m)- (9) 

Given a formalism to describe the observed CMB perturbations, it is then nec- 
essary to calculate a theoretical description of the perturbations. The starting 
point for this is to solve the Boltzmann equation for the radiation transfer 
of photons. To proceed, we expand the perturbations in Fourier modes of 
wavcvcctor k. A full derivation of the necessary equations is beyond the scope 
of this paper (for details see Refs. [21,32,43]), so we will summarise the im- 
portant equations below. 

Tensor perturbations are assumed to arise, in similar fashion to scalar pertur- 
bations, from quantum fiuctuations during infiation. Although our knowledge 
of this epoch is speculative, we may describe the perturbations by a primordial 
power spectrum, 

Ph{k) = Ark^^-', (10) 

where k is the comoving gravitational-wave wavenumber and At is an am- 
plitude fixed by the energy density during inflation. The form of this power 
spectrum is motivated by inflationary theories which generically predict the 
tensor spectral index nr ~ 0, a nearly scale-invariant spectrum. If the process 
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generating the perturbations is Gaussian, then this power spectrum encodes 
all information about the distribution. 



The evolution equations for the tensor modes may be derived from the Einstein 
equations and are [44] 

Xij + '^-Xij + k'^Xij = l&T^Ga^'Kij. (11) 
a 

Here, Xij is a symmetric, traceless perturbation to the spatial part of the met- 
ric, and TTjj is the tensor part of the anisotropic stress. Here and throughout, 
overdots denote derivatives with respect to conformal time, and a(T) is the 
scale factor normalised to unity today. 

The temperature and polarisation anisotropics induced by an equal mixture of 
tensor modes of + and x polarisation with amphtude h may be described in 
terms of the variables Ax(to, n, k), where X=(T,E,B). These have dependence 
on both (j) and 6*, which motivates the use of new variables AT{T,fi,k) and 
Ap{t,h, k) dependent only on 6, first introduced by Polnarev [32]. Here, //. ~ 
fi ■ k is the angle between the direction n of propagation of the photon and 
the wavevector k of the tensor mode. The relation between these two sets 
of variables is detailed in Ref. [38]. For our purposes it is enough to use the 
Polnarev variables for calculating sources; the sources are then simply related 
to the original variables. 

The evolution of a single Fourier mode k satisfies the Boltzmann equations. 

At + ikfj.AT ^ -h- k[AT -"if], (12) 

Ap + ikiiAp = -fi[Ap + (13) 



* = 



Ir Ir 3~ 3~ 6-r 
Atc] -\- — At9 ~\~ A^d — ~Apri + — Zapo — ApA 

10 ° 7 70 5 ° 7 70 



(14) 



Here, we have defined the differential cross section for Thomson scattering 
as k — aUeXeaT, where rig is the electron number density, is the ionisation 
fraction, and ar is the Thomson cross section. The total optical depth between 
a conformal time r and Tq is given by integrating k to obtain k{t, Tq) — 
k{T)dT. The multipole moments of temperature and of polarisation are 
defined by A{k,ii) = X);(2/ + 1)(— i)'A;(A;)Pi(/x), where P;(/u) is the Lcgendre 
polynomial of order I. This decomposition converts Eqs. (12) and (13) into an 
infinite hierarchy of equations connecting higher moments to lower moments. 
These equations have solutions [38] 



\Jl^.Jo drST{k,r)^, (15) 
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Aei = / dTSp{k,T) 

Jo 



-ji{x)+jl'{x) + 



+ 



X 



(16) 



/■TO 

As; = / dTSp{k,T) 

Jo 



X 



(17) 



where ji{x) is the spherical Bessel function. In these expressions, x — A;(ro — r). 
Defining the visibihty function g{T) — ke~'^, the sources are given by 



ST{k,T) - -he-'^ + g^, 



(18) 



and the power spectra by 

Cxx'i = i^nfj k^dkPh{k)Axi{k)Ax'i{k). 



(19) 



(20) 



It is straightforward to show that for statistically equal distributions of left 
and right circularly polarised gravitational waves, the TB and EE cross- 
correlations vanish. If there is a preference for either polarisation, then a 
non-zero TB and EE correlation will be observed [47]. 

Equations for the evolution of the ionisation fraction, and hence the optical 
depth, exist, but will not be dealt with here. Details of the relevant equa- 
tions and useful analytic approximations for the optical depth and ionisation 
fraction may be found in Ref. 



Where necessary we assume a fiducial ACDM cosmology with = 0.05, 
^DM — 0.25, Oa = 0.7, and with Hubble parameter parameterised hy h — 
0.72. 

The above set of equations forms the basis for our problem. Having written 
expressions for the power spectra, we must now exploit a mixture of physical 
and mathematical approximations to bring out their implications. 



3 A Tale of Tensor Modes 



Let us now try to obtain an intuitive understanding of how the features in the 
power spectra arise. This will help motivate the approximations that follow in 
later sections. Useful discussions of how polarisation is generated are given in 
Refs. [40,49,50]. 
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First let us discuss the temperature power spectrum. The temperature mul- 
tipole moments due to an individual gravitational-wave of wavenumber k ob- 
served at a conformal time Tq are 



with X = k{TQ — t). The second of the sourcing terms is localised to the 
surface of last scattering (SLS) by the visibility function; as a consequence of 
the restricted range this term is small and may be neglected at all angular 
scales. Between I — 200 and I — 800, the contribution from the g'^ term 
falls off more slowly than the integrated Sachs- Wolfe (ISW) term allowing it 
to become marginally relevant. At lower and higher /, the power generated 
by the second term dies off rapidly and is totally negligible (Fig. 1). The 
first term, which dominates this integral, involves an integral from the SLS 
to the present day. Its form tells us that the temperature power spectrum is 
sensitive to the evolution of gravitational waves from the SLS to today and not 
to the recombination history. This term is a form of integrated Sachs- Wolfe 
(ISW) effect that describes the anisotropy generated by changing gravitational 
potentials. We can understand this effect by recalling that a gravitational- wave 




Fig. L Comparison of power generated by the two source terms for temperature 

anisotropy. Plotted are the total power (solid line), ISW term only (dotted line), 
and term only (dashed line). The g'i> term is essentially negligible at all /. The 
normalisation here, and in all plots, is specified by setting At = l/8(47r)^ and 
riT = 0. 

alternately stretches and compresses space as it oscillates. A photon travelling 
past the gravitational-wave loses energy when its wavelength is stretched. 
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but gains energy when its wavelength is reduced. If the gravitational- wave 
amplitude evolves over the course of the oscillation, the photon will undergo a 
net change in energy. Tensor modes decrease steadily in amplitude and oscillate 
after horizon entry. As such, a photon travelling along the crest of a phase 
front will slowly gain energy as the overall amplitude of the gravitational- 
wave decreases. Photons travelling at an angle to the mode experience further 
red and blue shifting as they propagate through different phase regions. Their 
energy oscillates as a consequence. Between the SLS and today, the period- 
averaged amplitude of the tensor mode decreases, and so the mean energy of 
the photon increases. 

If we consider a late time, so that the amplitude of the tensor modes is es- 
sentially zero, then we see that the final energy of the photon is determined 
by whether it started its journey from the SLS at a trough or crest in the 
tensor mode. Photons starting at a crest will have gained more energy and 
appear hotter than average and vice versa for those starting at a trough. This 
simplistic picture is modified by the effect of power free-streaming from one 
angular scale to another as the Universe expands, which tends to smooth the 
resultant power spectra. 

The situation is very different for the polarisation anisotropies. These are 
generated by expressions of the form 



Here, the source is very firmly localised to the SLS and so is sensitive to 
the thermal history and gravitational-wave evolution at that time. This is 
sensible. Treating the early radiation bath as unpolarised (as we expect from 
suppression of anisotropy during the tightly-coupled regime) , then polarisation 
is generated by Thomson scattering of an anisotropic intensity distribution. 
Where does this anisotropy come from? In the rest frame of the scattering 
electron, photons arrive from all directions from a mean distance determined 
by the mean free path of photons at recombination. In propagating, these 
photons experience the ISW effect, discussed in the case of the temperature 
spectrum, and so arrive at the scatterer with altered temperatures. The result- 
ing anisotropic temperature distribution is scattered, generating polarisation 
which free-streams to the present epoch. 

In this way, we can understand the power spectrum. For modes with wave- 
lengths much larger than the horizon size at last scattering, incident photons 
experience very little ISW before the last scattering event and httle polarisa- 
tion is generated. Optimal ISW and thus maximal polarisation is generated 
by modes that enter the horizon at the time of penultimate scattering. The 
amplitude of the gravitational- wave decays most rapidly immediately on hori- 
zon entry (see Fig. 2) before setthng into oscillation with a slowly decreasing 




(22) 
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amplitude. Modes that enter the horizon before penultimate scattering lead 
to photons whose ISW samples this slowly decaying regime. Hence, they gen- 
erate significant polarisation, but less than for the optimal case. Note that the 
time between penultimate and last scattering will be about the width of the 
surface of last scattering. 
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Fig. 2. Evolution of a gravitational-wave. Wavenumber k satisfies kr^q = 10. Shown 
solutions are numerical without anisotropic stress (solid curve), numerical with 
anisotropic stress (dotted curve), radiation (long dashed curve), matter (short 
dashed curve) , and WKB (dot-dashed curve) . The two vertical lines denote t = 1/k 

and Teq. 

Translating this into the form of the polarisation power spectrum, we expect 
a slow increase in power at large scales peaking at the scale of the horizon 
at penultimate scattering. Immediately after this, we expect a large drop in 
polarisation corresponding to the transition between modes that enter the 
horizon between penultimate and last scattering and those that do not. Next, 
we expect a steady decline in power as modes have entered the horizon before 
penultimate scattering and so redshifted away before the ISW effect is gener- 
ated. This region will show a transition in slope between modes that entered 
the horizon in the matter- and radiation-dominated epochs. On scales smaller 
than the mean free path at recombination, the power will drop sharply as 
phase cancellation between differing crests and troughs becomes important. 

In this discussion, it is important to realise that only three scales have entered 
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the problem. These are the comoving horizon at recombination, the horizon 
at matter-radiation equahty, and the width of the last-scattering surface. Fig. 
3 shows how the features in the power spectrum correspond to these scales. 




10 100 1000 

1 

Fig. 3. Tensor power spectra. Curves from top to bottom are Cj^ , Cp^, and C^^. 
Vertical lines indicate important angular scales, from left to right: horizon at re- 
combination, Tji, horizon at matter-radiation equality, Tgg, and the width of the 
last-scattering surface, Atr. 



4 Gravitational-wave Evolution 

Expansion of the Universe leads to damping of the tensor modes as described 
by the term proportional to h in Eq. (11). This is the usual redshifting of radi- 
ation. In addition, the tensor modes may be sourced by anisotropic stress, Tiij. 
It has been shown [44,51] that anisotropic stress generated by free-streaming 
neutrinos acts to provide viscosity, further damping the tensor modes. This 
effect is important only while the energy density in neutrinos is a significant 
fraction of the total energy; i.e., during the radiation-dominated epoch. 

The tensor modes may be decomposed into two independent polarisation 
states, hx and h+. With this decomposition and a source term appropriate 
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for neutrino anisotropic stress, we have 



h, + 2-h, + k-'h, = -24/,(t) (441 f K[k{T - T')]h,{T')dT' , (23) 

where i = +, x, and = Pu/p with p the unperturbed density, and K{s) is 
given by 

, sins 3coss 3 sins , 

^(^) " - — + — • 

To a first approximation, wc may neglect the effect of anisotropic stress, though 
it should be included in detailed calculations. Without the source term, ana- 
lytic solutions for Eq. (23) in pure radiation and matter cosmologies may be 
expressed in terms of the spherical Bessel function j/(x), 

SI Tl nl' 

hUr) = h{0)Mkr) = /^(O)-^, (25) 



/imat(T)=3M0)^i^. (26) 

In a mixed radiation and matter dominated universe, the solution follows /irad 
initially before asymptotically becoming similar to /imat- The initial radiation 
dominated phase introduces a phase shift into /imat as now the boundary 
conditions do not preclude the spherical Neumann solution to the unsourced 
Eq. (23). When calculating the power spectra it is important to get this phase, 
which determines the peak positions, correct. This point was understood but 
not implemented in Ref. [45] and included implicitly by others [46]. 

The behaviour of these solutions is shown in Fig. 2 and splits into three main 
regimes. When kr 1, h evolves slowly and is approximately constant. Once 
/cT ~ 1, the amplitude decays away rapidly before entering an oscillatory phase 
with slowly decreasing amplitude, when /cr 3> 1. Physically, this corresponds 
to a mode that is frozen beyond the horizon until its wavelength is of order 
the horizon size at which point it enters the horizon and redshifts rapidly with 
the expansion of the Universe. 

Recombination occurs shortly after the Universe becomes matter dominated. 
For modes that enter the horizon during the matter-dominated regime and 
so have evolved little in the radiation-dominated epoch, we expect /imat to be 
a good description. For modes that entered during radiation domination, we 
expect that the transition from radiation to matter domination will affect the 
evolution significantly. 

The matter-radiation transition can be accounted for in a variety of ways. 
Most simple is to assume that the transition is instantaneous and to match 
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the amplitude and derivative of h on the boundary. This will be a good ap- 
proximation for waves with wavelength much longer than the time taken for 
the transition to take place. 



^instant * 



r < Tf 



{T^jT)[Aji{kT) + By^{kT)l T > Teq, 



(27) 



with 



A 



B = 



2 fv / eq 



i/cTeq COs(2/cTeq) + Sin(2/cTeq) 



eq 



2 — 2k^T^ — 2 cos(2A;req) — /cTeq sin(2A;req) 



(28) 



(29) 



Alternatively, we may consider the situation where the wavelength of the 
gravitational-wave is much shorter than the transition time. In this case, the 
gravitational- wave sees the background expansion vary slowly and a WKB ap- 
proach is appropriate. Ng and Spcliotopoulis [53] first presented this approach, 
although they were primarily interested in late time asymptotic limits and so 
neglected the behaviour near the classical turning point. Here we generalise 
their result making use of the uniform Langer solution for the WKB problem 
[52]. The result is 



h{T) 



T{kT) 



-1/4 



Soir, 



TV2(r + 2) V2 
X 1 2V7rC2Ai 



1/6 



|so(r 



2/3' 



+ VTrCiBi 



2/3' 



with 



and 



2s 



s + 2k 



(30) 



(31) 



ds 



Soir) = / ^T{s) -. (32) 

JkT ^ S 

Here, tt is the solution to r(A;r) = 0, Ai and Bi are Airy functions, and 
Ci and C2 are constant coefficients set by the boundary conditions h{0) = 
1 and h{0) = 0. For technical reasons, these boundary conditions must be 
extrapolated to small r via asymptotic approximation to Eq. (23) and then 
applied. Care must be taken in evaluating the above expressions when t > tt- 
These details are discussed further in Appendix B. This WKB expression 
reproduces the phase of h in both radiation and matter dominated regimes, 
but underestimates the amplitude. The close agreement between the WKB 
and anisotropic-stress curves in Fig. 2 is a numerical coincidence. 
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Other approaches exist to handle this transition from radiation to matter in 
a more pragmatic fashion [45,46]. 



We can get the scahng of h from a simple argument. Before horizon entry, 
the amplitude /i of a gravitational- wave is constant. After horizon entry, the 

gravitational-wave redshifts with the expansion as radiation and scales as h ^ 
1/a. Hence, the amplitudes of a gravitational- wave today and at horizon entry 
are related by /itoday/ Gentry = Ocntry/atoday Taking /icntry to bc independent of 
k, we have /itoday oc Oentry Horizou entry occurs when aentry-f^entry = k, and 
so from the scaling of H in the matter- and radiation-dominated epochs we 
obtain Ocntry oc when radiation dominated and Oentry oc when matter 
dominated. Thus, we obtain the scalings. 



h oc < 



1, k< 1/to, 

A;-^ 1/Teq > A; > 1/to, 

k ^ , A; > 1 / Tgq . 



(33) 



This result agrees with both the instantaneous-transition and WKB solutions 
when r 3> Teq. These scaling relations form the basis for scaling of the power 
spectrum. We expect l{l + l)Cf'^ to scale as [26,45] 



l(l + l)Cl' oc <^ 



1, I < Ir, 

i"'^, leq> i > Ir, 
I ^ k ^ ^eq* 



(34) 



It has been claimed [26,54] that there should be an extra region, I > Za, in 
which the width of the last-scattering surface becomes important and due to 
phase-damping the scaling goes as . However the dominant source of tem- 
perature anisotropy is the ISW effect, which is insensitive to the recombination 
history, and so we do not expect to see this behaviour in the temperature power 
spectrum. One way to see this is to examine the kernel in Eq. (15). On small 
scales, which enter the horizon before tr, the finite rise time of e~'^ alters 
the weight in the integral by a nearly constant factor. For all other modes, 
it is sufficient to simply truncate the range of the integral to between tr and 
To, effectively imposing instantaneous recombination. On the other hand, the 
polarisation anisotropy is generated near the SLS and will show phase can- 
cellation dependent on the width of the SLS. We will return to this point 
later. 



In the absence of reionisation, at low / the power spectra for the polarisation 
grow as [37]. Prom this and the above scaling arguments, we would expect 
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the power spectrum to scale as 



I < Ir, 



l{l + l)Cf ^ oc < 




(35) 



I > Ia, 



with X=(E,B). The effect of phase-damping extends to much lower / than 
would be indicated by these simple dimensional arguments. Consequently, the 
region of constant power l^q < I < Ia is never visible in calculated spectra, 
but is lost in the transition to the phase-damping regime. 

The above expressions for the gravitational-wave amplitude are used to gener- 
ate the power spectra displayed in Fig. 4. The plots are normalised by taking 
riT — and setting At = 1/8 (47r)^. All of the plots show the same scahng 
relation at low I. This regime is dominated by modes that have not entered 
the horizon at recombination and so are approximately constant. At low I 
hraat Underestimates the power, while /ij-ad overestimates the power. This is a 
consequence of the contribution of the modes that have entered the horizon 
that are evolving in a mixed radiation-matter universe and so have amplitudes 
intermediate to the predictions of these two approximations. 

Moving above the peak at we clearly see the different scaling relations be- 
tween the matter and radiation approximations. Recall that at these high Is, 
we expect the main contribution to come from modes that entered the horizon 
in the radiation-dominated epoch, and so h^ad should be a good approximation. 
The WKB result shows a transition between following the matter-dominated 
curve to behaving more like the radiation-dominated form, though with re- 
duced amplitude. This reduction in amplitude is an unfortunate characteristic 
of the WKB solution and is not significant to understanding the physics. 
The WKB solution serves as a nice bridge between matter- and radiation- 
dominated epochs. The instantaneous solution fails to be useful on scales with 
wavelength short compared to the transition time-scale. 

These curves display the scahng expected from Eq. (33), but we see that in 
the numerical case, excluding anisotropic stress, we never observe the full 
scaling for a matter-dominated regime. The combination of recombination oc- 
curring soon after matter-radiation equality and the Universe becoming mat- 
ter dominated only slowly means that the power spectra damp more slowly, 
closer to l"^, for scales 1/tr < k < 1/req. The presence of matter also causes 
peak positions to shift to smaller scales over the fully radiation-dominated 
case indicative of the phase shift that the transition introduces in h. 

For the purposes of reproducing the exact tensor-mode power spectra, we must 
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Fig. 4. T and B power spectra calculated using approximate forms for the gravita- 
tional-wave amplitude h. Plotted are the results using h from the full numerical cal- 
culation (solid curve) and from the radiation-dominated (long dashed curve), mat- 
ter-dominated (dot-short dashed curve), instantaneous-transition (dotted curve), 
and the WKB (dot-long dashed curve) approximations. 

worry about preserving both the amplitude and phase of the gravitational 
waves. The importance of the amplitude is clear in estimating the power cor- 
rectly. The phase determines the positions of the maxima and minima in the 
high-Z region of the spectra. Maxima correspond to gravitational waves whose 
amplitude was at a maximum or minimum at the SLS; minima correspond 
to gravitational waves whose amplitude was close to zero at the SLS. Alter- 
ing the phase of the gravitational waves shifts the k values for which these 
maxima occur at the SLS and so shift the features in the CMB. If we wish to 
understand these features in detail, then we must understand how the phase of 
the gravitational waves varies with k and how this is mapped onto the power 
spectrum. This mapping is the subject of Section 6. 



5 Recombination History 

While the Universe is young and hot, baryons are ionised and tightly coupled 
to photons via Thomson scattering. Once the temperature falls below a few 
eV, it becomes favourable for electrons and ions to recombine to form neutral 
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molecules. As the number of charged particles falls, the mean free path of any 
given photon increases. Eventually, the mean free path becomes comparable to 
the horizon size and the photon and baryon fluids are essentially decoupled. It 
is at this point in the Universe's evolution that the CMB photons last scatter. 

The visibility function describes the probability that a given CMB photon last 
scattered from a particular time. In terms of the optical depth k, this visibility 
function is given by 



Numerical calculations show that gir) is sharply peaked during recombination. 
This property suggests we approximate the visibility function by a narrow 
Gaussian for analytic simplicity. For example, 



determines the visibility function in terms of the conformal time tr of recom- 
bination, its width Atr, and the amplitude g{TR} at recombination. 

Approximating the visibility by a Gaussian leads to a simple form for the 
optical depth in the region close to tr. If we write n in the general form 
K — exp[— /(r)], then consistency with Eqs. (36) and (37) requires that k ~ 
exp[— (r — tr)/Atr] and g{TR) ^ l/{eATji) in the region close to tr. This 
latter result is essentially a statement about the normalisation of the Gaussian 
and preserves the total weight of the visibility function for different widths. 
Away from recombination, the evolution of the optical depth is a complicated 
function of the thermal history and not easily approximated. 

These approximations for g{r) and k are plotted in Fig. 5 for the fiducial 
cosmology with Qf, = 0.05, Qbm = 0.25, and Q\ = 0.7. For this cosmology, we 
have To = 13515 Mpc, tr/tq = 0.0203, t^Jtq = 0.0076, and Atr/tq = 0.0012. 
While the Gaussian form does a reasonable job of approximating the shape of 
the peak, the visibility function is clearly skewed and possesses a significant 
tail. The combination of these features means the Gaussian approximation 
will underestimate the power and shift features to slightly smaller angles than 
in the true power spectrum. 

Fig. 6 shows a series of power spectra calculated using the Gaussian approxi- 
mation. In each, the correct thermal history is used to calculate the evolution 
of the source function \E' with the Gaussian approximation applied when cal- 
culating the Ax from Eqs. (18) and (19). Although not strictly self-consistent, 
this isolates the modification of the source due to a changed thermal history 
from the effect of the visibility function on generating anisotropies. Source 
evolution will be considered in Section 7. 

The temperature power spectrum shows no variation with Atr at I < 200. 




(36) 



g{T) = g{TR)e 



(37) 
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Fig. 5. Recombination history. Plotted are the visibility function g{T) and the optical 
depth K calculated numerically (solid curves) and the approximations described in 
the text (dashed curves) using Atr = 15.7. 
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Fig. 6. Evaluation of the Gaussian approximation for the visibility function. Three 
Gaussian derived power spectra are shown for Atr = 10 (long dashed curve), 15.7 
(short dashed curve), and 22 (dotted curve). 
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Power on these scales is generated via the integrated Sachs- Wolfe effect by 
modes that only evolve significantly between and Tq and so arc insensitive 
to the thermal history. At smaller scales, the modes of interest are evolving 
over recombination and so contain information about the thermal history. 
Modifying the width of the visibility function affects the power spectrum via 
the term in Eq. (21) which acts to cut the integral off below tr. Widening 
the SLS makes this cutoff slower which, owing to the concave nature of e"'^, 
leads to less weight in the integral. This leads to the differences observed in 
the top panel of Fig. 6. This is not phase-damping, and does not alter the 
scaling of the power spectrum significantly. In addition to this overall shift in 
power, larger Ari? acts to wash out the bumps and wiggles. For modes that 
oscillate rapidly over this rise time, the ISW samples an averaged starting 
value of h and so is less sensitive to the presence of peaks and troughs. Power 
is still generated from the net decrease in the maximum amplitude of h from 
Tr to To, so this cffect does not cause a large drop in power as suggested in 
Refs. [26,54]. 

Fig. 6 shows that the Gaussian approximation leads to a lower polarisation 
power spectrum. This difference is a consequence of the long tail to the vis- 
ibility function, which is not reproduced in the Gaussian approximation. In 
these plots, varying Ar^ docs not affect the overall amplitude of the power 
spectrum. This is an artefact of using the same ^ for each plot. In reality, the 
amplitude of the polarisation power spectra depends sensitively on Ari^, as 
will be shown in Section 7. These plots show that varying Atr in the Gaussian 
approximation does not affect the shape at low Z, but a wider width leads to 
a sharper fall off in power at high I. This is a feature of phase-damping, which 
will be discussed in Section 7. 

None of the three values used precisely reproduces the dechne of the true power 
spectrum, which is seen to fall off more rapidly than the approximations. This 
seems to be a consequence of the tail to the visibility function. In keeping with 
expectation, the peaks in the high-/ region are found at slightly higher / in the 
approximations than the numerical result. 



6 Projection Factors 



The power spectra that we observe today are projections of the temperature 
and polarisation anisotropies at the last-scattering surface. By inspection of 
Eqs. (15), (16), and (17) we can define three projection terms, 

Pti{x) = (38) 
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Pei{x) = -Mx) + f/ix) + ^ + MM, (39) 

Pbi{x) = 2j[{x) + (40) 

Typically the argument of these terms is A;(ro — r), the look-back time scaled 
by the wavenumber, reflecting that these are projections from the point of 
origin onto today's sky. 

The different forms of the projection factors, plotted in Fig. 7, help explain 
many of the features seen in the power spectra (Fig. 8). 
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Fig. 7. Projection terms for I = 50. Dashed curve shows approximations. 





First consider the E projection factor as a simple example. A sharp peak occurs 
at X ^ I. This tells us that the value of Cf^ at / = 50 is determined by the 
behaviour of the source function at a; ~ /. The polarisation source function is 
strongly peaked around r = tr, which implies the behaviour at last scattering 
of the mode with wavenumber k ~ 1/{tq — tr) dominates the contribution 
to Ci. If the projection factor was a Dirac delta function, this would be the 
whole story. However, the projection factor has a significant tail for x > I 
signifying that modes with larger wavenumber also contribute power to this 
angular scale. 

From this, we can see that the sharper the spike at x ^ I, the sharper the 
features seen in the power spectrum. A wider peak mixes in modes of different 
phases blurring the spectra. Noting that the B projection factor lacks a sharp 
peak, we expect the B power spectra to contain blurred features relative to 



19 




Fig. 8. Power spectra generated solely due to the source term g^. Temperature 
(solid curve), E mode (dotted curve), and B mode (dashed curve). 

the E spectra. In addition, as its maximum is at higher x, we expect features 
in the source to be shifted to smaller / than in the E spectrum. This sort of 
argument has some validity with the T spectrum, but is there complicated by 
the extended nature of the source term. 

The complicated form of the projection factors makes analj^ic progress dif- 
ficult. Similarly, their oscillatory behaviour makes numerical integration very 
slow at high k values. CMBFAST implements a scheme for fast numerical 
integration. Here we discuss time averaging the projection functions to get a 
useful analytic envelope. 

Widely known approximations for the spherical Bessel functions in the cases 
X 3> / and X <^ I exist and are in common usage. For our purposes, though, 
we are most interested in the case where a; ~ /; i.e., in the vicinity of the peaks 
of the projection factor. A relatively simple approximate form may be derived 
which is valid in the regime x > I [55], 

ji{x) — . : COS [x(sinQ; — a cos a) — 7r/4] , (41) 

V x^ sin a 

where cos a = {I + l/2)/x. This can be shown to reduce to the usual ji{x) ~ 
sin(x — lTr/2)/x for x ^ Approximations valid in the regime x ~ I exist, 
but are more complicated and will not significantly improve on this level of 
approximation. 

This approximation is still complicated and shows strong oscillation. When we 
calculate power spectra we will be interested in quantities of the form [Px{x)]^. 
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To proceed, we substitute Eq. (41) into the projection factor and then average 
the squared projection function over a full cycle to extract the variation of 
the envelope. Time averaging makes use of the relations (sin^ x) = (cos^ x) — 
1/2 and (sin a; cos x) = 0. This envelope may then be further simplified by 
assuming x ^ 1 and / ^ 1. The resulting envelope functions are not pretty, 
but may be numerically integrated by a standard routine. They are 



{PE{xf) ~ { - 16(/ + f- x^){l + 12/=* + 8/^ + Sx" - A{x + 2lxff 
+ (-16/^(2 + /) - 4(1 + 2/)(3 + m)x^ + 32a;^ 
+ (-1 + 8/(1 + 2l)){x + 2lxff}/{bl2x\-l{l + /) + x^'^}, (43) 

{12f-Sx^f~m + 2f-2x^ni + f-x^) 
^^^^^^ ^ x3(-4F + 4x2)5/2 ■ ^44^ 

The form for {Pt^xY) is consistent with that quoted in Ref. [56]. 

One slight complication is that these approximations show divergent behaviour 
as X — * / making it necessary to arbitrarily restrict the domain to ,x > / + a, 
where a is an arbitrary cutoff of order unity. Fig. 7 shows these approxima- 
tions and the cutoff. This cutoff procedure behaves best in the case of the B 
projection factor which is already decreasing as x — > L The E and T projection 
factors have considerable weight near to the peak, and so care must be taken 
in selecting the cutoff. 

To apply these projection factors, we approximate the anisotropy term. 

Ax; = r dTg{T)'^{T)Pxi[k{T^ - r)], (45) 

by the expression, 

Axi ~ PxiHto - tr)] r rfr^(r)^(r), (46) 





where we have pulled the projection factor out from the integral. This should 
be a good approximation provided that the projection factor varies slowly 
relative to the source term. 

The power spectra calculated under this approximation are shown in Fig. 9. 

Note that two approximations are combined here. The projection term has 
been pulled from the integral and then the time averaged expression used in 
place of the exact form. The inadequacy of this approximation at large / can 
be seen. This is to be expected, as in this regime the relevant k modes oscillate 
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Fig. 9. Comparison of power spectra with exact (solid curve) and approximate 
(dashed curve) projection factors. 

significantly over the width of the last-scattering surface. At values I < 500, 
the approximation is appropriate. 

A significant discrepancy in the EE power spectrum is visible at / < 6. In 
making the approximations to produce Eqs. (42), (43), and (44), we have 
discarded information about very large angular scales. The rise in power in 
the EE power spectrum is sensitive to these large angular scales [38] and so is 
not reproduced by our approximations. 



7 Source Evolution 

Having discussed the recombination history and the mechanics of projection, 
we turn to the core issue of how the source itself evolves. 

Those modes that enter the horizon close to recombination evolve only slowly 
over the width of the last-scattering surface. For these modes, it is possible 
to derive analytic approximations to the source functions that occur in the 
expressions for Ax- Here, we follow the approach of Zaldarriaga and Harari 
[56], but see also Ref. [57]. Preceding recombination, the optical depth n is 
large and the photons are tightly coupled to the baryonic fluid. In this regime, 
we may expand the Boltzmann equations for the temperature and polarisation 
multipoles in powers of k~^. Keeping terms to first order in k,~^, we obtain 
the equations, 

Aro = -/i-«[Aro-*], (47) 
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Apo = -k[Apo + % 



(48) 



Atz = , Z > 1, (49) 

A« = , / > 1. (50) 

Using these equations together with the definition of ^ gives us an expression 
for the time evolution of the source function within this tightly-coupled limit, 



Thus, an approximate solution for ^ is 

h{r') 



K{r') — k{t) 



10 



(52) 



If we assume that the visibility function is approximately Gaussian during 
recombination, then we may approximate k ~ —k/Atr. This allows a change 
of variable to x — k{t')/k{t) which leads to 

*(r) = -^e^'^i^^Arn f ^e'^-. (53) 

In taking the gravitational-wave driving term outside of the integral, we have 
assumed that h varies slowly over the visibility function. This approximation 
is only valid for k ^ 1/Atr. At larger wavenumbers, the rapid oscillation of 
\l/ over the visibility function makes this a poor approximation. An improve- 
ment is to replace h{TR) in Eq. (53) with its value averaged over the visibility 
function 

{h{T)) = r dTg{T)h{T) ^ /i(r«)e-(^^-«)V^ (54) 
Jo 

In calculating the right-hand side, we have treated h as an oscillatory func- 
tion with a slowly-varying envelope. Integrating an oscillatory function over 
a Gaussian leads to the function evaluated at the Gaussian's peak multiplied 
by a decaying exponential. This exponential decay has a clear physical inter- 
pretation. For modes with k > I/Atr, the source function oscillates rapidly 
across the visibility function. Hence, different regions in the SLS contribute 
to the observed polarisation with different phases leading to cancellation and 
a decrease in the observed power. We will refer to this cancellation as phase- 
damping. While present in the scalar modes, this effect is overwhelmed by 
diffusion damping [58]. Diffusion damping makes the effective visibility func- 
tion scale dependent and always sufficiently narrow that phase-damping is 
not important. For the tensor modes, phase-damping provides the dominant 
process for damping on small scales. 
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Fig. 10. Comparison of Eq. (53) (long dashed curve), Eq. (53) with damping term 
(short dashed curve), and numerical calculation of ^' (solid curve). Four different 
values of k are plotted: Utr = 0.03, 1.23, 12.2, and 28.2. An arbitrarily scaled 
visibility function (dotted curve) has been plotted in each panel to guide the eye. 

Fig. 10 shows the behaviour of ^ and our analytic approximations for four 

values of k. For small k, the approximation closely mirrors the growth of 
in the region where the visibility function has weight. At larger k, the source 
function is seen to oscillate across the width of 5'(t); this is not reproduced by 
either approximation. This should not be cause for concern as we now discuss. 

The quantity of real interest is the anisotropy that this source generates. This 
is calculated from Eqs. (16) and (17). For the given polarisation, X=(E,B), 
we have 

/■TO 

Axi(k)^ drg(r)^(r)Pxi[k(ro-r)]. (55) 
Jo 

First we pull the projection term outside of the integral assuming that it varies 
slowly over the width of the visibility function, 

/•TO 

Axi{k)^Pxi[k{To-TR)] dTg{T)<i/{T). (56) 

Jo 

Note that ^ appears only through an integral over the visibility function. 
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k (Mpc-i) 

Fig. 11. Comparison of J^' for the numerical * (solid curve), the approximation 
for ^ without phase-damping (dotted curve) and the approximation for ^ with 
phase-damping (dashed curve). Vertical lines indicate (from left to right) I/tr, 
1/rcq, and I/Atr. 

Provided that our approximation can reproduce this integrated behaviour, 
the fact that it fails to reproduce the temporal oscillation is unimportant. Fig. 
11 indicates the close agreement between the integrated source function and 
our approximation, provided that phase-damping is taken into account. 

Having checked the validity of our approximation we substitute for ^ in (56) 
using Eqs. (53) and (54) giving 

^xi{k) = PxiUto - T«)]-^/i(r«)Ar^e-('=^-«)'/2 H dKe'^o- T ^e"^-. 

lU Jo Ji X 

(57) 

The integrals evaluate to (10/7) log(10/7), which leads to the final result 

Ax; = Pxi[k{T^ - t«)]A(t«) AT«e-(^^^«)'/^ Q log . (58) 

This result is proportional to the width ^tr of recombination as might be 
expected. During recombination, photons will travel for a distance of order 
Atr before scattering. This is the time available for the quadrupole which 
sources the polarisation to grow, and so we expect a result proportional to 

Extending this result to calculate the power spectrum is straightforward. We 
have Eq. (20), 

Cxi = (47r)2 j k^dkPhik) [Axi{k)f . (59) 



25 



Applying our expression for Axi{k) yields the final result for this Section, 
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Fig. 12. Comparison of E and B power spectra for CMBFAST (solid curve) and 
analytic approximations: h for radiation epoch (short dashed curve), h calculated 
numerically (dotted curve) , and h calculated numerically with exponential damping 
(long dashed curve). 

Fig. 12 compares the result of CMBFAST with that from Eq. (60). AU plots 
have been calculated for the fiducial cosmology with = 0.05, JIdm = 0.25, 
and = 0.7. Values for and At/j were chosen by fitting by eye to the 
visibility function produced by CMBFAST. Radiation and numerical forms for 
h without phase-damping are plotted alongside a numerical form for h with 
phase-damping and the results of the full numerical calculation. 



Agreement between all of the solutions is good at low / where our assumptions 
are most valid and the anisotropy is building slowly. Similarly, the position of 
the main peak is accurately reproduced, although it becomes clear that the E 
projection-factor approximation is less reliable than the B-mode one. Beyond 
the main peak, the effects of the different forms for h become apparent. 

Without phase- damping, the radiation-dominated form of h leads to an almost 
fiat power spectrum, consistent with Eq. (35), with a pattern of bumps and 
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wiggles roughly in phase with those in the full numerical calculation. Moving 
to the numerical form of h leads to some damping over the radiation case. This 
reflects the increased redshifting which occurs when matter becomes important 
between Teq and tr. We never observe the expected matter-dominated scaling 
of although the presence of matter does lead to a shght dechne in power as 
/ increases. For this cosmology, the ratio p,-/ Pm evaluated at tr is 0.29 showing 
that the radiation content is still significant at recombination. 

Neither of these forms reproduces the rapid decline in power at large /, which is 
not surprising as we have yet to include phase-damping. Once this is included, 
the shape of the power spectrum is much closer to that of the numerical 
calculation showing a sharp decline in power above / ^ 300 and reproducing 
the position of the peaks to reasonable accuracy. Power in the range / = 150 to 
/ = 600 is slightly overestimated. For / > 600, the limitations of our projection- 
factor approximations become apparent with the very rapid drop in power 
previously observed in Section 6. 



8 Discussion 



Here, we discuss the information that the features in the tensor power spec- 
trum contain and how detection would complement our existing understanding 
of the early Universe. In contrast to the scalar modes, the tensor modes con- 
tain very clean information about the evolution of the Universe. The features 
of the scalar spectra are a result of the oscillation of the matter-radiation fluid 
during the period up to recombination. The scalar spectra encode informa- 
tion about the sound speed of the baryon-radiation fluid, the baryon fraction, 
and other information about the particle content of the Universe [49,59]. In 
contrast, the features of the tensor spectrum are determined solely by the 
wave motion of the evolving gravitational waves. They primarily contain in- 
formation about the expansion rate during the early Universe. Through their 
overall angular scale, both spectra encode basic information about the epoch 
and duration of recombination and the geometry of the Universe. 

The first peak of the tensor polarisation spectrum, occurring at / ^ 90 for the 
B modes and / ~ 105 for the E modes, is determined by the horizon scale at 
recombination. The exact angular scale is determined by this along with the 
redshift of recombination and the geometry of the Universe. This information 
can be determined from the scalar modes allowing a direct measurement of the 
horizon scale. The amplitude of the tensor power spectra is directly related to 
the energy scale of inflation. Slow-roll inflation, parameterised by the energy 
scale Ei of inflation, predicts a B-mode power spectrum with a peak at Z ~ 90 
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and peak amplitude [60] 

A5pea. = 0.024 (^^A_)%K. (61) 

Measuring this main peak is the subject of several experimental endeavours 
[27,28,30] with hope of detection in the not-so-distant future. 

For all models consistent with WMAP constraints to the energy scale of in- 
flation, features after the main peak are sub-dominant to the lensed B-mode 
signal [60,61]. This necessitates the use of algorithms to clean the polarisation 
maps and recover the tensor signal. Techniques using maximum likelihood 
[62,63,64] and quadratic estimators [65,66] have been advanced to deal with 
this problem. Even so, this will complicate precision measurements of the ten- 
sor B mode after the main peak. 

Measuring the overall amplitude after the main peak should recover the scaling 
relations discussed in Section 4. The breaks in the different regimes yield the 
horizon scales of matter-radiation equality and the width of the SLS. This 
in itself is enough information to constrain a cosmological model yielding 
and fir- Neutrino anisotropic stress further damps power on small scales and 
will make detection more difficult while adding extra information about the 
neutrino fraction. 

The positions of the peaks and troughs on small scales contain information 
about the phase of the gravitational-wave at recombination. This in turn de- 
pends upon the early expansion rate. The acoustic peaks in the scalar power 
spectrum can be used as a standard ruler; the wiggles in the tensor power spec- 
trum can have the same utility, but operating on a different range and spacing. 
Measuring these wiggles would better constrain our cosmology, though it is 
doubtful that scales small enough to probe the very early Universe will be 
observed in the foreseeable future. 



9 Conclusions 

This investigation has probed the individual elements that compose the cal- 
culation of the tensor power spectra. Using a variety of approximations, we 
have obtained a semi-analytic expression which qualitatively reproduces the 
behaviour of more detailed calculation. While the approximation is clearly not 
suitable for precise comparison with data, it serves to illustrate the important 
physics in an intuitive fashion. We have shown that the features of the power 
spectrum may be explained with reference to three main scales: the horizon 
size at recombination, the horizon size at matter-radiation equality, and the 
width of the SLS. The first two scales determine the evolution of the tensor 
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modes; the latter relates to the effect of the thermal history on the generation 
of anisotropies. The shapes of the polarisation spectra show most sensitivity 
to the width of the SLS through the phenomenon of phase-damping that dom- 
inates on small scales. The effect of the thermal history on the temperature 
spectra is much less dramatic, affecting the amplitude and smoothing on small 
scales. We have seen that the position of the peaks and troughs in the power 
spectrum relate to the phase of the gravitational- wave at recombination. It is 
of interest that we do not see modes displaying matter- dominated behaviour 
in our calculations. This would allow the tensor spectra to probe the radiation 
content at recombination. Useful scaling relations have been developed and 
clarified. We hope this paper will aid in general understanding of the tensor 
modes and inspire future experimental efforts. 



A Numerical Evolution of Gravitational Waves vi^ith Anisotropic 
Stress 

In this Appendix, we return to the question of anisotropic stress. In the early 
Universe, free-streaming neutrinos provide the main source of anisotropic 
stress. After recombination, photons free stream and can also contribute, 
though the energy density in radiation is falling fast and the effect is negligi- 
ble. Working from Eq. (11) and standard expressions for the energy density of 
a distribution of relativistic massless particles, an integro-differential equation 
describing the evolution of the tensor modes may be derived [51], 

h + 2-h + eh = -24/,(r) ( r K[k(T - T')]h(T')dT', (A.l) 

where fi, = p^/p with p the unperturbed density, and K{s) is given by 

, sins 3coss 3 sins ,. 

The new term acts to damp the amplitude of h and can be seen to have the 
form of a convolution over the mode's past history of the kernel K{s) and 
the "velocity" h. The hnear dependence on means that the damping term 
will become negligible in the matter-dominated regime where fi, oc a~^. In the 
radiation-dominated epoch, though, the neutrino and total energy densities 
scale in the same way leading to f^, = 0.40523. This suggests that the damping 
term will primarily affect those modes that enter the horizon well within the 
radiation-dominated epoch. In consequence, it will affect the power spectrum 
only at high I where these modes are the dominant contributors. 

The right-hand side is also damped by the {a/a}'^ term which scales as in 
both the matter- and radiation-dominated epochs. 
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Numerical integration of Eq. (A.l) is possible after recasting the integro- 
differcntial equation as a set of coupled Voltcrra integral equations [67]. These 
coupled equations may then be integrated using standard techniques [68]. 

Given an integro-differential equation of the form, 

y^'^Xt) ^f{t,y{t),..., y^'-'\t)) + 1^ K(t, s, y{s), y^^\s))ds, (A.3) 

and defining Zk{t) = y^''\t) for A; = 0, . . . , r — 1, we may recast Eq. (A.3) as 
the set of first-order Volterra integral equations, 

zr{t) = J^K (t, s, y{s), . . . , y^'-'\s)) ds, (A.4) 

Zr-l{t) = "^^ + (/(S, Zo{s), Zr-l{s)) + Zr{s))ds, (A.5) 

Zk{t)=yo'^+[ Zk+i{s)ds, (A.6) 
Jo 



where the y^ are the relevant initial conditions. 

Making the co-ordinate transformation t — kr in Eq. (A.l), we can apply this 
formalism to obtain 

Z2{t) - ['Knit - s)zi{s)ds, (A.7) 
Jo 

zi{t) = yl + [-2^z^{s) - zo{s)) + Z2{s)^ ds, (A.8) 

Mt) ^yo+ f zi{s)ds, (A.9) 
Jo 

where primes denote differentiation with respect to t, and we have defined 
Zi — h', Zq — h, and 

K^{s) = -2AMt) (^^) K{t - s). (A.IO) 



Fig. 2 shows the results of numerical integration of Eq. (A.l). The inclusion 
of anisotropic stress damps the wave during the radiation-dominated regime 
leading to decreased amplitude and a slightly shifted phase. 

Of primary interest here is the effect of the anisotropic stress by the time of 
recombination. It is traditional to calculate a transfer function relating the am- 
plitude and phase of the numerical solution to that of the matter-dominated 
solution hjaa.t{T) — 3ji(/cr)/A;T. At recombination, radiation is still important 
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and this analytic form is a relatively poor approximation. For illustrative pur- 
poses, we will numerically calculate the amplitude ratio A and phase shift \1/ 
between numerical calculations of h with and without the effects of anisotropic 
stress evaluated at the time of recombination. To calculate A at a given r, we 
first calculate and then numerically fit /istress(^'?') with ^4 /ino-stress(^'?' + ^) 
over the period containing r. This avoids the oscillation that results if we seek 
to obtain the amplitude ratio by simply dividing /igtressl^T) by /ino-stressl^T). 

The results of this calculation (Fig. A.l) illustrate that anisotropic stress intro- 
duces a /c-dependent damping asymptoting to a factor of A ~ 0.81. The phase 
shift introduced remains small and reaches a maximum value of \E' ~ 0.13 rad. 
The amplitude ratio is unity while the gravitational-wave remains within the 

1.05 

1 

0.95 
^ 0.9 
0.85 
0.8 

0.15 
0.1 
0.05 


Fig. A.l. Transfer function for gravitational waves when anisotropic stress is present. 
Top panel shows the variation of the amplitude ratio A between the case with 
anisotropic stress and without anisotropic stress. Bottom panel shows how the phase 
shift ^ between the two cases varies with k. 

horizon. Around horizon entry its rises slightly above unity before decreasing 
asymptotically to A ~ 0.81. This slight rise is a consequence of the anisotropic 
stress which, like viscosity acting on a pendulum, slows the initial decrease of 
the gravitational- wave and leads to a lower final amplitude of oscillation. The 
phase difference between damped and undamped cases grows after horizon 
entry as a consequence of the slower evolution of the damped wave. It peaks 
and begins to asymptote to a constant value for modes that entered the hori- 
zon sufficiently before matter-radiation equality to reach their asymptotically 
damped form. 

CMBFAST may be modified to incorporate the evolution equation (A.l) and 
the resulting power spectra calculated. Fig. A. 2 shows the effects for the T and 




31 



B power spectra. The spectra are essentially unchanged at low /. Only at high 
/ do we see a suppression of power from the additional damping. Damping only 
occurs in waves that have evolved significantly during the radiation-dominated 
epoch, so this makes sense. 




10 100 1000 

1 

Fig. A. 2. T and B power spectra incorporating anisotropic stress (dashed curve) 
and without anisotropic stress (solid curve). Damping in the power is clearly seen 
on small scales. 



B WKB Solution 

If the scale factor changes more slowly than the evolution of the gravitational- 
wave, then WKB techniques become a sensible method of approximation. In 
this Appendix, we detail the application of this approach to the evolution 
of a gravitational-wave mode through the matter-radiation transition. The 
WKB approach was first applied in this context by Ng and Speliotopoulos 
[53] , although the solution presented here is of a slightly more general nature. 

We begin with the equation of motion for /i(r) and specific initial conditions 
that we wish to solve. Working with the dimensionless variables rj = (\/2 — 
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l)'''/'^eq s-nd q — kTeq/ {^2 — 1), we have 



h + 2-h + q^h^O, 
a 



(B.l) 



h{0) = 1, 



(B.2) 



h{0) = 0, (B.3) 

with overdots indicating differentiation with respect to rj. The behaviour of 
the scale factor in a universe containing only dust and radiation is given by 



i{v) = aeciViv + 2) 



(B.4) 



where is the scale factor at equality. 



To move towards the standard WKB form, we make the transformation h 
y/[viv + 2)] which eliminates the first derivative term leading to 



y+{q- 



0. 



(B.5) 



Before proceeding, we notice that this transformation is singular aX rj = 0, the 
point at which our boundary conditions arc specified. This prevents us from 
applying the boundary conditions (b.c.) directly and will motivate looking for 
asymptotic solutions to (B.l) which we will discuss later. 

Wc next map from the interval [0, oo] to [—00, 00] by the transformations 
Tj = cxp(x') and u{x) = exp{—x/2)y{x). These place our problem in the WKB 
form 

u"ix)^f{x)u{x), (B.6) 



The standard WKB problem [52] 

eV{x) = Q{x)y{x), 



(B.7) 



(B.8) 



where Q{0) — 0, Q{x) > for x < and Q{x) < for x > 0, has the uniform 
Langer solution 



/ 3 \ 



X <^ 2CoAi 



3 . ^ 



+ CiBi 



3 . ^ 



(B.9) 
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with 



'S'o = /° jQit)d-t- 

Jx 



(B.IO) 



Here, Ai(x) and Bi(x) are Airy functions and Ci and C2 are constant coeffi- 
cients to be set by tfie boundary conditions. 

Applying this directly to Eqs. (B.6) and (B.7) and manipulating the algebra 
slightly, we obtain 



/iwkb(T) = —- -So{t) 



TV2(r + 2) V2 
X J 2C2Ai 



o ^ 2/3- 

5S„(r) 



with 



and 



r(.) 



2s 



1 

4 ^ s + 2k 



S,{t) ^ / yrXs) -. (B.13) 

JkT * S 

Here, tt is the solution to Vikr) — 0. It can be shown that Utt is bounded 

such that 1/2 < krT < 3/2. These expressions may be evaluated directly 
when r < tt. but some care must be taken when r > tt- In this case, r(A;r) < 
0, and we must keep careful track of minus signs. Making the appropriate 
manipulations, we obtain for t > tt 



+ CiBi 



o X 2/3- 



fB.Ul 



(B.12) 



ds 



Vi[-r(fcr)]-V^ /3 
rV2(r + 2) [2'^^' 



with 



and 



X <^ 2C9Ai 



2/3' 



+ CiBi 



2/3' 



2s 



4 S + 2A; 



ds 
s 



, (B.14) 

(B.15) 
(B.16) 



Asymptotic expansions of the Airy functions recover the more familiar expo- 
nential and trigonometric forms for the WKB connection formula. 

Having obtained an expression for h, we now need to apply the boundary 
conditions. Unfortunately, if we naively try to apply the boundary conditions, 
we discover that /iwkb('?" = 0) is divergent. This is a consequence of the trans- 
formation required to place the equation of motion into WKB form which is 
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singular at r = 0. To get around this problem, we seek to apply the bound- 
ary conditions at some small time when the gravitational- wave has had little 
chance to evolve and will be well described by a series solution. 

Eq. (B.l) has a regular singular point at r = so we attempt a Probenius 
series solution of the form 



(B.17) 



n=0 



Substituting this expression into Eq. (B.l), solving the indicial equation for 
cr, and equating like powers of r to get a recurrence relation leads to an 
asymptotic polynomial expression for h. In fact, this procedure only generates 
one of the two linearly independent solutions to Eq. (B.l). Application of the 
boundary conditions causes the other solution to vanish and normalises this 
one leaving us with a solution. 



h{r) = 1 - + 



6 
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240 



(3 - 2q^)T\ 



(B.18) 



This is valid only when r ^ 1. We can use this to extrapolate the b.c. from 
zero time to some small time and use it to determine the constant coefficients 
in Eq. (B.U). 
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Fig. B.l. Comparison of the WKB (dashed curve) and numerical (solid curve) evo- 
lution of h for k = 0.1677Mpc^^ in the fiducial cosmology. The amplitudes /imat 
(long dashed curve) and /irad (dotted curve) are plotted for reference, and Teq is 
indicated by a vertical line. 

The WKB solution (Fig. B.l) accurately reproduces the phase of the oscillation 
in the regime r > tt although it underestimates the amplitude of h in this 
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region by a factor of ~ 0.87. The approximation is also good for t < tt failing 
only at times comparable to the time at which the b.c. are applied. For the 
power spectra calculated in this paper, this is taken to be /cTmatch = 10~^. 
Consequently, the power spectrum become unreliable at small multipoles, / < 
10. On these scales power is generated by modes that have evolved very little 
by the time of recombination. 
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